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Only  a  couple  months  have  elapsed  since  the  last  annual  report,  so  there  is  nothing  significant 
to  report  for  those  intervening  months.  The  technical  accomplishments  and  outcomes  from  the 
grant  support  have  been  documented  in  the  three  annual  reports,  which  are  attached  as  appendices. 
We  provide  here  only  summary  statistics  for  the  three  year  period,  and  a  summary  of  the  most 
significant  accomplishments. 

Noteworthy  accomplishments  include 

•  development  of  genetic  algorithms  (GAs)  involving  both  discrete  and  continuous  variables  for 
composite  structure  design; 

•  Fortran  90  templates  and  fully  distributed  massively  parallel  implementations  for  genetic  al¬ 
gorithms; 

•  validation  of  the  superiority  of  the  normal  flow  homotopy  algorithm  for  tracking  nonlinear 
equilibrium  curves; 

•  use  of  response  surface  approximations  to  optimal  surfaces  in  global/local  optimal  design; 

•  use  of  compatibility  metrics  at  the  global  level  to  improve  blending  between  optimized  local 
subsystems; 

•  binary  tree  memory  and  local  approximations  with  continuous  variables  to  improve  GA  effi¬ 
ciency; 

•  parallel  GAs  with  migration  for  blending  panels  in  the  design  of  large  composite  aircraft 
structures; 

•  cellular  automaton  approaches  to  structural  design; 

•  decomposition  theory  for  quasi-separable  multidisciplinary  design  optimization  problems. 

The  grant  supported  the  publication  of  8  journal  papers  and  28  conference  papers.  During 
the  grant  period  an  additional  12  journal  papers  were  accepted  and  28  more  journal  papers  were 
submitted  for  publication.  Four  graduate,  four  undergraduate,  and  six  AASERT  students  were 
supported,  with  one  MS  ^ilesis  completed  and  three  Ph.D.  theses  in  progress.  Three  Pis— Zafer 
Giirdal,  Raphael  T.  HaM^,  Layne  T.  Watson — ^were  partially  supported  by  the  grant.  The  project 
hosted  five  visiting  scholars  and  involved  three  post-doctoral  researchers.  There  were  five  distinct 
(different  contacts)  industrial  technology  transfers.  The  PI  Watson  received  the  Alumni  Award  for 
Excellence  in  Research  in  2000. 
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Objectives 


Global-local  design  of  large  structures  with  complex  geometry  presents  a  great  challenge  to 
structural  designers.  In  particular,  coordination  between  codes  that  optimize  the  structure  or 
major  structural  components  (such  as  entire  wing  or  fuselage  structures)  and  codes  that  optimize 
individual  panel  details  requires  substantial  and  careful  effort  to  reach  optimal  or  near  optimal 
solutions.  Moreover,  the  compatibility  of  the  adjacent  designs  (or  the  lack  of  it)  in  terms  of  their 
geometric  and  material  variables  presents  a  serious  manufacturing  difficulty  for  the  large 
structure.  We  are  presently  developing  a  computational  infrastructure  and  algorithms  with  a 
sound  theoretical  basis  to  extend  industrial  ad  hoc  approaches  to  the  global-local  design  and 
blending  of  local  designs.  We  proposed  a  two-level  optimization  approach  employing  genetic 
algorithms  tailored  to  panel  design  on  the  lower  level.  Response  surface  approximations  to 
optimized  panel  failure  loads  are  then  used  for  the  upper  level  wing  or  fuselage  optimization.  In 
addition,  solutions  to  the  blending  problem  at  the  local  level  using  genetic  algorithms  will  be 
investigated.  We  also  proposed  the  development  of  a  novel  kind  of  population  based  cellular 
automata,  similar  to  genetic  algorithms,  for  enforcing  construction  compatibility  between 
adjacent  panels.  Parallel  computing  and  high-level  object-oriented  code  development  will  be 
employed  at  every  opportunity. 


Status  of  Research 

This  document  summarizes  the  work  performed  from  September  1,  1998  to  August  30,  1999. 
During  this  period,  progress  was  made  on  the  development  of  a  global/local  design  capability  for 
structures  exhibiting  geometrically  nonlinear  behavior  and  on  refining  capabilities  for  wing 
design  under  strength  and  buckling  constraints. 

The  work  on  geometrically  nonlinear  behavior  was  motivated  by  engineers  at  the  Boeing 
Company,  who  were  interested  in  optimizing  fuselage  structures  for  postbuckling  response. 
Traditional  design  procedures  assume  that  a  structure  fails  once  buckling  Cccurs.  In  reality,  many 
structures,  such  as  plates  and  stiffened  panels,  are  able  to  safely  carry  a  significant  amount  of 
additional  load  beyond  the  buckling  load  (see  Figure  1).  Structures  that  are  capable  of  safely 
carrying  this  additional  load  are  said  to  have  postbuckling  strength.  If  the  structure  is  designed  by 
taking  into  account  this  postbuckling  strength,  significant  weight  savings  can  be  obtained 
compared  to  traditional  buckling  based  designs.  In  weight  critical  applications,  such  as  fuselage 
design,  designers  often  attempt  to  utilize  the  structure’s  postbuckling  strength  to  save  weight. 
Unfortunately,  no  affordable  optimization  procedure  incorporating  postbuckling  analysis 
currently  exists  to  assist  these  designers  in  their  task.  Existing  procedures  are  either 
computationally  prohibitive  or  inaccurate.  The  goal  of  this  effort  is  to  use  the  global/local 
methodologies  previously  developed  under  this  project  to  develop  and  demonstrate  an  affordable 
design  methodology  for  designing  structures  for  local  postbuckling  response. 


First,  the  capability  to  perform  affordable  design  optimization  of  postbuckled  composite  plates 
and  stiffened  panels  was  developed.  This  was  achieved  by  connecting  a  computationally  efficient 
geometrically  nonlinear  analysis  code  to  a  genetic  algorithm  computer  code.  Improvements  in 


the  computational  efficiency  of  the  design  process  were  obtained  by  implementing  the  normal 
flow  algorithm  for  solving  the  nonlinear  equations  describing  the  response  of  the  structure. 
Preliminary  global/local  results  were  generated  for  the  design  of  a  simple  composite  fuselage 
model. 


Figure  1  -  Postbuckling  Response  of  Blade  Stiffened  Panel 


The  second  thrust  in  the  project  was  to  refine  the  coordination  process  between  global  and  local 
design  in  the  linear  range.  This  was  applied  to  wing  structural  design,  also  motivated  by  interest 
from  Boeing  (formerly  McDonald  Douglas)  engineers  and  researchers.  One  focus  of  the  wort 
was  to  improve  the  efficiency  of  the  process  of  obtaining  a  large  number  of  panel  designs  via 
genetic  algorithms.  This  was  achieved  through  introducing  repair  operators  into  the  genetic 
algorithm  used  for  panel  design.  The  repair  operators  deal  with  constraints  much  more  efficiently 
than  the  traditional  treatment  via  penalty  techniques. 

A  second  focus  in  coordination  was  to  improve  the  accuracy  of  the  response  surface  of  optimal 
panel  failure  loads  used  for  communications  between  the  global  and  local  design  optimization.  A 
normalization  procedure  was  implemented  and  currently  we  are  exploring  the  use  of  derivatives 
of  the  panel  failure  loads  for  improving  further  the  accuracy.  Derivatives  have  not  been 
traditionally  used  in  response  surface  construction,  because  response  surfaces  were  originally 
developed  for  fitting  experimentally  obtained  data.  Consequently,  very  little  theory  is  available 
on  efficient  response  surface  construction  using  derivative  data.  One  present  effort  is  to  develop 
the  required  theoretical  basis. 

Finally,  students  supported  on  the  AASERT  part  of  the  research  sought  the  application  of  the 
same  local-global  design  philosophy  to  the  design  of  Micro  Aerial  Vehicles  (MAV).  Previous 


work  under  the  project  sought  the  coordination  of  wing  and  airfoil  design.  However,  the  coupling 
between  these  two  aspects  of  the  aerodynamic  designs  proved  to  be  too  weak  to  demonstrate  the 
benefits  of  the  approach.  Current  work  is  focused  on  coordination  of  the  overall  wing  design  and 
the  design  of  the  control  system  needed  to  provide  the  MAV  with  sufficient  roll  stability.  MAVs 
can  carry  very  little  weight,  and  the  design  of  the  control  system  can  substantially  add  to  the 
required  battery  weight  if  the  vehicle  is  difficult  to  control.  Therefore,  the  design  of  the  shape  of 
the  vehicle  is  closely  coupled  with  the  design  of  the  control  system. 


Global/Local  Design  Formulation 

The  global/local  design  procedure  for  postbuckling  response  will  be  first  demonstrated  for  a 
fuselage  design  example.  At  the  global  level,  a  finite  element  model  of  a  fuselage  segment  is 
used  (see  Figure  2).  The  goal  of  the  global  optimizer  is  to  minimize  the  fuselage  weight  without 
allowing  structural  failure  to  occur.  Global  design  variables  include  the  thickness  of  0°,  ±45°,  and 
90°  plies,  and  constraints  are  imposed  on  the  stresses  in  the  skin  and  frame  elements.  The 
individual  fuselage  skin  segments  are  each  designed  using  a  local  postbuckling  analysis.  Local 
design  variables  control  the  stacking  sequence  of  the  skin,  and  constraints  are  imposed  on  the 
collapse  load  and  material  failure. 


Figure  2  -  Global  Fuselage  Model 


Design  Capability  for  Local  Postbuckling  Response 

Before  the  design  of  the  overall  structure  begins,  the  local  optimizer  is  used  to  construct  a 
database  of  optimized  local  designs.  Each  design  in  the  database  is  generated  as  a  function  of  the 
global  design  variables  and  the  applied  loads.  For  a  given  set  of  global  design  variables  and 
applied  loads,  the  local  stacking  sequence  will  be  designed  to  maximize  the  failure  load.  Since 
the  global  design  variables  and  the  internal  loads  vary  as  a  function  of  the  global  optimization 
iterations,  the  local  design  database  will  be  constructed  for  a  range  of  these  values.  The  range 


will  correspond  to  realistic  combination  of  those  global  variables  dnd  loads  that  may  be  observed 
in  a  realistic  structure. 

Once  the  local  design  database  is  constructed,  the  design  of  the  overall  structure  can  begin.  For 
any  set  of  global  design  variables  and  applied  loads,  the  global  optimizer  can  consult  the  local 
design  database  for  each  skin  panel  to  determine  the  failure  load  of  the  optimized  panel.  In  this 
way,  the  global  optimizer  obtains  information  about  the  optimized  local  designs  and  can 
incorporate  this  information  into  the  design  of  the  global  structure  without  actually  doing  the 
local  optimization  in  real  time. 

In  order  to  implement  a  global/local  design  methodology  for  postbuckled  designs,  the  capability 
to  optimize  local  panels  in  the  nonlinear  regime  must  be  developed.  This  can  be  accomplished  by 
linking  a  computationally  efficient  nonlinear  analysis  code  with  an  appropriate  optimization 
code.  In  this  project,  a  nonlinear  finite  strip  analysis  was  adapted  to  communicate  with  a  genetic 
algorithm  based  optimization  code.  The  finite  strip  analysis  [1],  which  was  developed  as  part  of  a 
previous  research  project,  is  capable  of  approximating  the  geometrically  nonlinear  response  of 
compressively  loaded  prismatic  plate  structures  such  as  plates  and  stiffened  panels.  The  code 
includes  the  effects  of  geometric  shape  imperfections  and  is  capable  of  identifying  and  traversing 
limit  points  in  the  response.  The  genetic  algorithm  [2],  which  was  developed  earlier  under  the 
present  grant,  allows  us  to  design  the  stacking  sequences  of  plate  and  stiffener  elements  while 
imposing  realistic  manufacturing  constraints. 

During  the  development  of  the  nonlinear  design  capability,  the  design  process  was  made  as 
computationally  efficient  as  possible.  To  this  end,  a  new  algorithm  for  solving  the  nonlinear 
system  equations  was  implemented  in  the  finite  strip  analysis.  This  algorithm,  the  normal  flow 
algorithm,  was  found  to  yield  significant  computational  savings  compared  to  the  previously 
implemented  solution  algorithm  (a  variant  of  the  Crisfield  method).  For  several  of  the  test 
problems,  the  normal  flow  algorithm  required  fewer  than  half  the  number  of  matrix 
factorizations  in  order  to  tra^e  the  nonlinear  response  up  to  the  point  at  which  the  structure 
collapsed.  In  an  effort  to  call  more  attention  to  the  normal  flow  algorithm  in  the  structural 
mechanics  community,  a  iechnical  paper  was  written  and  submitted  to  the  International  Journal 
of  Solids  and  Structures  [3]. 

In  order  to  increase  further  the  computational  efficiency  of  the  nonlinear  design  procedure,  the 
combined  analysis/design  code  was  modified  so  that  multiple  local  optimization  runs  could  be 
performed  in  parallel.  This  capability,  which  is  still  being  tested,  will  be  of  great  benefit  in  the 
generation  of  the  data  required  for  building  the  local  response  surfaces. 


Design  Study 

In  order  to  demonstrate  the  global/local  methodology  and  the  potential  benefits  of  designing  for 
postbuckling  response,  a  simple  design  study  will  be  performed.  The  fuselage  segment  illustrated 
in  Figure  2  will  be  designed  for  several  different  value  of  the  maximum  postbuckling  load  factor 
(PFmax).  The  postbuckling  load  factor  (PF)  is  defined  as  follows: 

PF=1-P,/P, 


where  Pt  is  the  skin  buckling  load  and  Pf  is  the  skin  failure  load.  The  maximum  postbuckling 
load  factor,  PFmax.  is  the  maximum  value  of  the  postbuckling  load  factor  that  is  allowed  to  occur 
for  a  given  design.  This  factor  governs  how  much  of  the  total  applied  load  is  allowed  to  be 
carried  in  the  postbuckling  regime.  For  example,  if  PFmax  is  0.33,  the  skin  panels  will  only  be 
allowed  to  carry  33%  of  the  total  load  in  the  postbuckling  regime.  If  PF  is  0.0,  none  of  the 
applied  is  allowed  to  be  carried  in  the  postbuckling  regime  (this  corresponds  to  a  traditional 
buckling  based  design).  In  this  project,  fuselage  designs  will  be  generated  corresponding  to 
several  different  values  of  PFmzx  in  order  to  demonstrate  the  potential  advantages  of  designing 
for  postbuckling  response. 

Preliminary  results  corresponding  to  PFmax  —  0.33  were  generated  for  several  portions  of  a 
fuselage  segment  obtained  from  Boeing.  For  the  purposes  of  the  design  study  described  here, 
however,  a  simplified  fuselage  model  is  being  used  (Figure  2).  This  simplified  model  will  allow 
the  results  of  the  design  study  to  be  published  in  the  open  literature,  and  will  allow  other 
researchers  to  duplicate  our  results.  This  simplified  fuselage  model  is  assumed  to  be  loaded  by 
bending  loads  only,  resulting  in  a  simple  biaxial  stress  state  for  each  skin  panel.  Locally 
optimized  designs  for  the  simplified  model  are  currently  being  generated.  Initial  designs  will 
neglect  the  effects  of  load  redistribution  at  the  global  level  due  to  postbuckling  deformations.  If 
these  effects  are  significant,  they  will  be  aceounted  for  later  in  this  project. 


Global-Local  Wing  Design  Coordination 

The  improved  global-local  design  coordination  was  investigated  in  the  context  of  global-local 
wing  design.  The  overall  wing  modeled  was  analyzed  and  designed  using  the  GENESIS  finite 
element  based  structural  optimization  code.  GENESIS  is  a  commercial  code,  in  use  in  both  the 
automotive  and  aerospace  industries.  It  employs  continuous  design  variables,  so  one  of  the 
hurdles  that  we  needed  to  address  was  the  process'  of  rounding  the  continuous  designs  obtained 
from  GENESIS  for  use  in  the  panel  design  optir  uzation  that  employs  discrete  variables.  It  was 
found  that  simple  rounding,  followed  by  min^:  adjustments  in  the  design  were  sufficient  to 
address  this  problem. 

The  global-local  design  procedure  includes  the  following  stages:  (i)  generation  of  large  number 
of  local  panel  designs;  (ii)  fitting  a  response  surface  to  the  failure  loads  of  the  optimal  designs; 
(iii)  overall  wing  design  employing  the  response  surface;  and  (iv)  panel  design  using  loads  and 
overall  thicknesses  obtained  from  the  overall  wing  design.  The  procedure  was  successfully 
demonstrated  on  a  simple  wing  stmcture  with  up  to  54  independent  global  design  variables. 

The  next  step,  now  in  progress  calls  for  integrating  derivative  information  in  the  creation  of  the 
response  surface.  Derivatives  of  the  failure  loads  of  the  optimal  panel  with  respect  to  the  global 
design  variables  and  loads  are  quite  inexpensive  to  obtain.  Therefore,  calculating  these 
derivatives  and  using  them  in  the  creation  of  the  response  surface  can  be  used  to  improve 
accuracy  at  low  cost  or  to  reduce  computational  cost  at  comparable  accuracy.  The  improved 
accuracy  was  demonstrated  derivatives  of  panel  loads,  and  work  is  in  progress  on  obtaining  and 
using  derivatives  with  respect  to  global  thickness  variables. 


Parallel  GAs  for  Blended  Local  Optimal  Designs 


In  previous  work  a  Fortran  90  GA  framework  was  developed  for  use  with  composite  laminate 
design  optimization.  The  framework  included  a  GA  module,  encapsulating  GA  data  structures 
and  basic  operations,  and  a  package  of  GA  operators,  which  worked  with  the  module.  An 
extension  to  the  framework  was  then  introduced.  The  distributed  GA  extension  implemented  a 
migration  algorithm  to  work  with  a  group  of  sub-populations  evolving  in  parallel.  The  goal  of  the 
migration  algorithm  was  to  improve  the  normalized  cost  and  reliaMIiiy  of  a  set  of  GA 
optimization  runs  without  significantly  impacting  the  time  it  takes  to  make  those  runs.  We  are 
building  on  this  research  in  a  new  direction  now,  considering  the  benefits  of  communication 
between  populations  of  different  optimizations  running  in  parallel  potentially  to  improve  the 
blending  of  the  local  structural  components. 

An  MPI  based  master-slave  Fortran  90  program  was  developed  to  allow  us  to  collect  response 
surface  data  for  the  local/global  optimization  in  parallel.  This  simple  program  is  the  first  step  in 
developing  a  far  more  useful  tool  to  study  how  to  promote  continuity  of  properties  between  GA 
optimized  adjacent  panels.  The  global  optimization  provides  constraints  on  the  lower-level 
design  allowing  the  local  GA’s  to  operate  on  a  fixed  number  of  plies  of  each  orientation.  The 
problem  then  becomes  reconstructing  the  stacking  sequence  for  each  local  optimization  as  a 
permutation  of  these  plies.  Currently,  each  GA  run  is  conducted  locally  with  no  knowledge  of 
the  structure  of  the  surrounding  panels.  Local  designs  created  in  this  way  often  have  little  in 
common  with  adjacent  panels  causing  manufacturing  costs  to  rise.  Fortunately  GA's  can  develop 
many  near  optimal  solutions  with  widely  different  stacking  sequences  for  the  same  local  design. 
We  are  exploring  methods  for  these  local  optimizations  to  run  in  parallel,  communicating 
through  structure  inherited  topologies  and  migration,  to  construct  a  design  that  is  optimized  for 
both  weight  and  manufacturing  cost. 
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Objectives. 

Global-local  design  of  wing  and  fuselage  structures  presents  a  great  challenge  to  structural 
designers.  In  particular,  coordination  between  codes  that  optimize  the  entire  wing  or  fuselage 
structure  and  codes  that  optimize  individual  panel  details  requires  substantial  and  careful  effort 
to  reach  optimal  or  near  optimal  solutions.  We  are  presently  working  with  a  group  at  Sikorsky 
Aircraft  and  another  group  at  Boeing  Long  Beach  (formerly  McDonnell  Douglas)  on  solutions  to 
this  problem.  This  proposal  calls  for  developing  a  computational  infrastructure  and  algorithms 
with  sound  theoretical  basis  to  extend  industrial  ad  hoc  approaches.  We  propose  a  two-level 
optimization  approach  employing  genetic  algorithms  tailored  to  panel  design  on  the  lower  level. 
Response  surface  approximations  to  optimized  panel  failure  loads  are  then  used  for  the  upper  level 
wing  or  fuselage  optimization.  We  also  propose  the  development  of  a  novel  kind  of  population 
based  cellular  automata,  similar  to  genetic  algorithms,  for  enforcing  construction  compatibility 
between  adjacent  panels.  Parallel  computing  and  high-level  object-oriented  code  development  will 
be  employed  at  every  opportunity. 


Status  of  research. 

Global/local  optimization  and  blending  using  genetic  algorithms  with  migration: 
Aircraft  stiffened  composite  panels  are  typically  designed  by  optimizing  them  locally,  resulting 
in  incompatibilities  between  the  stacking  sequence  of  the  adjacent  panels  that  form  the  entire 
structure.  Such  incompatibilities  often  result  in  structures  that  are  difficult  to  manufacture,  or 
in  some  cases  impossible  to  manufacture.  Although  the  differences  are  driven  by  local  loads  and 
constraints  that  are  applied  across  the  structure  at  different  panels  differently,  it  may  be  possible 
to  improve  the  manufacturability  by  directly,  or  indirectly,  enforcing  the  adjacent  laminates  to 
have  similar  stacking  sequences.  In  particular,  when  genetic  algorithms  are  used  to  optimize  the 
local  panel  stacking  sequences,  often  multiple  near  optimal  solutions  with  minor  differences  in  the 
laminate  stacking  sequences  are  obtained  for  each  panel.  Therefore,  it  may  be  possible  to  accept 
small  losses  from  optimality  of  the  adjacent  panels  but  achieve  substantial  gains  in  the  continuity 
of  the  laminate  stacking  sequences  between  adjacent  panels  may  be  possible.  In  implementing  the 
genetic  algorithm  optimization,  by  allowing  local  designs  to  evolve  in  parallel,  sending  migrants 
adjacent  panels,  the  objective  functions  may  be  intelligently  modified  to  favor  a  globally  blended 
design.  During  the  present  effort  the  blending  of  the  laminate  stacking  sequence  is  accomplished 
using  the  string  edit  distance  between  individuals  and  the  set  of  emigrants  from  adjacent  panels. 
Parallel  results  are  presented  for  multiple  communication  topologies  as  well  as  varying  migration 
and  blending  rates. 

In  the  present  research  design  of  composite  material  aircraft  structures,  particularly  designing 
the  local  details  of  stiffened  panels  in  the  wing  and  fuselage.  These  panels  are  made  of  composite 
laminates,  using  fiber-reinforced  layers.  The  design  of  the  fiber-reinforced  composite  laminate 
requires  the  specification  of  the  stacking  sequence  with  orientation  and  material  type,  creating 
a  discrete  optimization  problem.  It  is  computationally  expensive  to  design  an  entire  wing  or 
fuselage  structure  with  the  panels  optimized  simultaneously.  Instead,  local  panels  are  commonly 
optimized  (failure  load  is  maximized  given  a  fixed  number  of  plies  of  each  orientation)  for  the 
specified  local  loads.  An  improvement  to  this  process,  a  global/local  procedure,  fits  a  response 
surface  to  the  maximal  failure  load  of  the  local  panel  configuration  as  a  function  of  the  local  loads 
and  the  number  of  plies  of  each  orientation.  The  response  surface  is  then  used  in  a  global  wing 
optimization  (with  numbers  of  plies  of  each  orientation  as  design  variables)  providing  constraints  on 


2 


panel  load  carrying  capacity  while  taking  into  account  the  variation  of  the  local  load  redistribution 
as  a  function  of  the  design  changes.  The  global  solution  (the  load  distribution)  is  then  used  one 
final  time  to  construct  locally  optimal  panels  that  match  the  overall  wing  design.  The  stacking 
sequences  of  the  final  panel  designs  in  this  scenario  are  usually  incompatible  across  the  overall 
structure. 

The  global/local  optimization  strategy  briefly  described  above  yields  designs  that  are  locally 
optimal  but  potentially  infeasible  as  a  final  solution  due  to  manufacturing  incompatibilities  between 
adjacent  panels.  The  incompatibilities  themselves  result  from  the  large  differences  in  stacking  se¬ 
quence  (the  difference  in  the  fiber  orientation  angle  at  different  through-the-thickness  locations  in 
the  laminate)  across  the  overall  structure.  The  research  conducted  during  the  present  reporting 
period  resulted  in  the  development  of  an  algorithm  that  takes  advantage  of  the  parallel  implemen¬ 
tation  of  the  genetic  algorithm  optimization  to  enforce  some  level  of  continuity  in  the  laminate 
stacking  sequence  between  adjacent  panels. 

The  advantage  of  the  algorithm  proposed  here  lies  in  the  nature  of  multiple  elitist  selection 
in  genetic  algorithms  themselves.  In  a  standard  elitist  selection  strategy  only  a  single  member  of 
a  parent  population  can  survive  the  selection  process  and  be  placed  in  the  child  population.  In  a 
multiple  elitist  selection  strategy  the  genetic  algorithm  allows  for  a  greater  number  of  high  fitness 
members  to  survive  the  selection  process  at  each  generation.  Building  on  the  fact  that  a  local  panel 
optimum  (or  near  optimum)  can  appear  in  several  distinct  stacking  sequences,  perhaps  a  globally 
blended  solution  can  be  derived  via  intelligent  selection  from  the  elite  of  a  population. 

The  particular  focus  of  this  work  is  in  using  genetic  algorithms  in  parallel  to  achieve  compatible 
local  panels.  Each  individual  panel  is  designed  by  a  different  processor  in  a  parallel  processing  envi¬ 
ronment.  The  populations  representing  each  panel  evolve  in  parallel  and  send  migrant  individuals 
to  adjacent  populations.  The  migrant  population  carries  information  about  the  best  individuals 
from  adjacent  populations  to  compute  string  edit  distances  from  the  working  population  to  the 
migrant  population.  The  computed  distance  value  is  used  to  reward  those  individuals  in  a  popula¬ 
tion  that  are  evolving  closely  to  those  found  in  the  migrant  populations  from  adjacent  panels.  This 
re  .vSirding  process  is  a  small  modification  to  the  objective  function  fitness  value  and  serves  to  blend 
panel  designs  of  the  population  balancing  the  optimization  of  weight  with  stacking  sequence 
compatibility.  Comparative  results  are  shown  for  a  base  line  of  oblivious  unblended  panels  with 
those  panels  obtained  from  multiple  migration  and  blending  variations  as  well  as  the  effect  of  varied 
communication  topologies  on  the  propagation  of  information  and  communication  time. 

Geometrically  nonlinear  analysis  of  local  panels:  During  local  panel  optimization,  it  is 
necessary  to  compute  the  nonlinear  relationship  between  the  applied  loads  and  the  resulting  struc¬ 
tural  deformations.  By  tracking  this  nonlinear  equilibrium  path,  it  is  possible  to  identify  critical 
buckling  and  collapse  phenomena.  In  structural  mechanics,  the  nonlinear  equilibrium  equations 
are  often  solved  using  the  Newton  Raphson  method  in  conjunction  with  an  incremental/iterative 
solution  method  (Haisler  and  Stricklin,  1972).  In  order  to  trace  the  equilibrium  paths  through 
limit  points  in  the  response,  the  “arc  length”  methods  of  Riks/Wempner  (Wempner,  1971  and 
Riks,  1979)  or  Crisfield  (Crisfield,  1981)  are  commonly  employed.  A  closely  related  method,  which 
is  based  on  a  “normal  flow”  algorithm  (Watson  et  al.,  1987  and  Watson  et  al.,  1997),  is  less  popular 
despite  the  fact  that  it  may  have  advantages  in  terms  of  computational  efficiency  and  robustness. 
During  the  reporting  period,  an  implementation  of  the  normal  flow  algorithm  to  nonlinear  lo¬ 
cal  composite  panel  analyses  is  investigated,  and  the  potential  advantages  of  the  method  for  the 
solution  of  structural  mechanics  problems  is  demonstrated. 
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In  the  following,  a  brief  description  of  the  Riks/Wempner,  Crisfield,  and  normal  flow  solution 
methods  is  provided.  Numerical  results  are  presented  for  the  collapse  of  a  blade  stiffened  panel.  The 
relative  efficiencies  of  each  method  are  compared  and  the  advantages  of  the  normal  flow  algorithm 
in  this  case  are  illustrated. 

The  discrete  set  of  nonlinear  equilibrium  equations  representing  the  response  of  a  structural 
system  may  be  written  in  the  following  general  form: 

f(A,d)  =  0,  (1) 

where  f  €  is  a  nonlinear  function  of  a  scalar  loading  parameter  A  and  the  displacement  vector 
d  e  In  an  incremental/iterative  solution  method,  these  equations  are  solved  in  a  series  of 

steps  or  increments,  usually  starting  from  the  unloaded  state  (A  =  0) .  The  n-th  step  begins  from 
a  known  solution  on  the  equilibrium  path,  =  (A;,d;)^,  and  consists  of  a  prediction  phase  and 
a  correction  phase.  In  the  prediction  phase,  an  estimate  for  the  next  point  on  the  equilibrium 
path,  =  (A°+i,d°_f.i)^,  is  generated.  Beginning  at  this  point,  Newton  Raphson  iterations 
are  employed  during  the  correction  phase  to  find  a  new  point  on  the  equilibrium  curve.  The 
Riks/Wempner,  Crisfield,  and  normal  flow  algorithms  discussed  here  all  utilize  an  approximation 
to  the  arc  length  of  the  equilibrium  path,  ,  to  fix  the  size  of  each  step.  The  i-th  iteration  executed 
by  each  algorithm  at  step  (n  +  1)  can  be  described  using  the  same  general  set  of  equations: 


Dt  (zrA)]Az  =  -f(4+i). 

(2) 

^ (^n+l)  ~  0) 

(3) 

where 

Az  =  z^^.1  —  i  =  l,2,.... 

and 

r*  (Zn+ijj  dX  dd 

is  the  m  x  {m  +  1)  Jacobian  matrix  of  the  system  (2)  .  Equation  (3)  is  an  auxiliary  “constraint” 
equation  that  is  different  for  each  algorithm.  Upon  convergence  of  the  Newton  Raphson  iterations, 
a  new  solution,  z*+i,  is  obtained  and  the  prediction/ correction  process  is  continued  for  subsequent 
steps. 

It  should  be  noted  that  the  system  of  equations  (2)-(3)  are  not  usually  solved  simultaneously, 
as  equation  (3)  destroys  the  symmetric  and  banded  qualities  that  equation  (2)  usually  possesses. 
Instead,  a  procedure  similar  to  that  first  proposed  by  Batoz  and  Dhatt  (1979)  is  often  employed 

(Ramm,  1980).  If  the  modified  Newton  Raphson  algorithm  is  used,  Di  (z^+\)  and  f  (zn+\)  are 
updated  only  at  the  beginning  of  each  step. 

Riks/Wempner  algorithm:  In  the  Riks/Wempner  algorithm,  the  auxiliary  equation  (3)  takes 
the  following  form: 

c(4„)  =  Kf  (<+!-<) -A!  =  0,  (4) 

where  z*  is  the  unit  tangent  to  the  converged  solution  at  step  n.  This  equation  defines  a  hyperplane 
that  is  normal  to  z*  and  which  is  at  a  “distance”  Al  from  the  previously  obtained  solution  at  step 
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n.  During  the  Newton  Raphson  iterations,  the  successive  iterates  are  forced  to  return  to  the 
equilibrium  path  along  this  hyperplane. 

CrisSeld  algorithm:  In  the  Crisfield  algorithm,  the  auxiliary  equation  defines  a  hypersphere  of 
radius  Al  centered  on  the  converged  solution  at  step  n: 

^  (^n+l)  =  (^n+1  “  ^n)  (^n+1  ~  ^n)  ~  ~ 

Successive  iterates  are  confined  to  the  surface  of  this  hypersphere  as  they  converge  to  the  equilib¬ 
rium  solution. 

In  practice,  A  is  sometimes  omitted  from  equation  (5),  and  the  constraint  is  instead  imposed 
in  m-dimensional  space  (Crisfield,  1981): 

c  (4+i)  =  (4«  -  <)’■  (4+1  -  4)  -  A'"  =  0-  («) 

This  confines  the  iterates  to  lie  on  a  cylinder  in  (A,d)  space  parallel  to  the  A  axis. 

Normal  flow  algorithm:  In  the  normal  flow  algorithm,  successive  Newton  Raphson  iterates 
converge  to  the  equilibrium  solution  along  a  path  which  is  normal  (in  an  asymptotic  sense)  to  the 
so-called  Davidenko  flow.  The  Davidenko  flow  can  be  described  by  considering  a  small  perturbation 
S  to  the  nonlinear  system  equations: 

f(A,d)  =  5.  (7) 

As  the  perturbation  parameter  varies,  small  changes  will  occur  in  the  solution  curve  for  equation 
(7).  The  family  of  curves  generated  by  varying  6  is  known  as  the  Davidenko  flow  (Allgower  and 
Georg,  1990). 

The  Newton  Raphson  iterate  Az  for  the  normal  flow  algorithm  is  the  unique  minimum  norm 
solution  of  the  m  x  (m  +  1)  equations  (2).  This  solution  may  be  obtained  in  two  steps.  First, 
a  particular  solution  v  to  the  equations  can  be  f"  imd  by  selecting  an  auxiliary  equation  (3)  an^i 
solving  the  resulting  system  of  equations  (2)-(3).  As  long  as  this  system  of  equations  has  rank 
m  -h  1,  the  auxiliary  equation  may  be  chosen  arbitrarily.  Once  a  particular  solution  is  obtained, 
the  minimum  norm  solution  Az  is 

v^u 

Az  =  V - ^u, 

u 

where  u  is  any  vector  in  the  kernel  of  [Df].  A  convenient  choice  for  u  in  this  case  is  the  unit 
tangent  vector  z* ,  which  has  usually  already  been  calculated. 

Numerical  results:  In  order  to  quantitatively  compare  the  relative  performances  of  the  solution 
algorithms,  a  sample  problem  was  solved  using  each  method.  This  sample  problem  consisted  of 
computing  the  response  of  a  compressively  loaded  graphite-epoxy  blade  stiffened  panel  loaded  into 
the  postbuckling  regime.  A  “stiffener-unit”  representation  of  the  panel  was  utilized,  wherein  the 
complete  panel  was  modeled  using  a  single  repeating  element  of  the  structure  comprised  of  a  single 
stiffener  and  the  adjacent  skin  (Figure  1).  The  stiffener  unit  was  48.1in  long,  6.28m  wide,  and 
had  a  stiffener  height  of  1.10m.  The  skin  layup  was  (±45/03)s,  the  stiffener  layup  was  (90/0)^, 
and  the  ply  thickness  was  O.OOGin.  The  skin  laminate  was  assumed  to  have  an  imperfection  in  the 
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Figure  1.  StifFener-unit  representation  of  panel. 

shape  of  seven  half-waves  in  the  longitudinal  direction  with  a  maximum  imperfection  amplitude  of 
0.0001m. 

Uniform  end  shortening  was  imposed  at  the  longitudinal  edges  of  the  structure,  and  symmetry 
boundary  conditions  were  imposed  at  the  lateral  edges.  At  the  longitudtiial  edges,  the'  skin  was 
simply  supported  out  of  plane  but  the  stiffener  was  not  allowed  to  rot'=t-  in  its  own  plane.  The 
panel  was  analyzed  using  a  finite  strip  analysis  of  the  semi-analytical,  multi^term  type  (Dawe,  1995 
and  Ragon,  1998).  The  finite  strips  were  modeled  as  balanced  and  symmetric  laminated  composite 
materials  which  were  assumed  to  behave  orthotropically  in  bending.  Material  properties  were  as 
follows:  El  =  20.0  •  lO^lb/in'^,  E2  =  1.30  •  lO^lb/in'^,  Gn  =  1.03  •  lO^lb/in^,  and  u  =  0.30. 

All  three  of  the  solution  algorithms  were  implemented  using  the  full  Newton  Raphson  method. 
The  prediction  phase  of  the  Riks/Wempner  and  Crisfield  algorithms  (including  the  determination 
of  the  arc  length  step  size  at  each  increment)  were  implemented  as  described  in  Crisfield,  1981 
under  the  heading  “A  Modified  Riks  Method” .  This  prediction  scheme  is  briefly  described  here. 

The  first  step  is  the  selection  of  an  appropriate  value  for  the  arc  length,  AL  The  arc  length  is 
adjusted  from  one  step  to  the  next  using  the  following  simple  formula; 

m 

where  m  is  the  number  of  iterations  that  were  required  at  the  previous  step  and  fh  is  the  (user 
specified)  desired  number  of  iterations  at  each  step.  This  procedure  allows  larger  steps  to  be  taken 
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Figure  2.  Normal  flow  trace  of  equilibrium  curve. 

when  the  solution  is  converging  easily,  and  forces  the  solver  to  take  smaller  steps  when  convergence 
is  more  difficult.  For  the  present  work,  fh  was  selected  between  3  and  5  so  as  to  achieve  the  most 
favorable  results.  Once  a  value  for  the  arc  length  has  been  computed  at  step  n  +  1,  the  following 
equation  is  used  to  predict  AA^^j: 

I  f 

where  Ad°_|,i  is  taken  to  be  the  tangential  solution  from  the  previously  converged  point. 

The  prediction  strategy  used  for  the  normal  flow  algorithm  was  that  implemented  in  subroutine 
STEPNX  in  HOMPACK90.  This  strategy  utilizes  a  Hermite  cubic  interpolating  polynomial  through 
previously  converged  points  and  a  sophisticated  step  size  estimation  algorithm.  The  load/end- 
shortening  behavior  of  the  structure  as  obtained  using  the  normal  flow  algorithm  is  illustrated 
in  Figure  2,  where  is  plotted  against  a  scalar  loading  parameter,  A.  A  value  of  A  =  1.0 
corresponds  to  a  uniform  axial  strain  =  —2.34  •  10“^  and  a  uniform  transverse  strain  in  the  skin 
of  =  +1.40  •  10~^.  The  response  is  linear  up  xmtil  Nx  ~  1130  lb,  where  the  skin  buckles  into  seven 
half-waves  along  the  length  of  the  stiffener-unit.  Beyond  this  point,  the  panel  continues  to  carry 
load  up  to  a  maximum  of  Nx  «  1620  lb,  where  the  panel  collapses.  The  normal  flow  algorithm  was 
able  to  continue  tracking  the  equilibrium  path  through  the  collapse  point  with  no  difficulty;  the 
algorithm  was  eventually  terminated  at  the  user’s  request  at  the  point  shown  (the  algorithm  did 
not  terminate  because  of  convergence  difficulties). 
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Figure  3.  Riks/Wempner  trace  of  equilibrium  curve. 

The  Riks/Wempner  and  Crisfield  algorithms  computed  the  same  load/end-shortening  curve 
as  did  the  normal  flow  algorithm,  but  they  both  encountered  difficulties  in  tracing  the  complete 
curve.  In  Figure  3,  a  zoomed  in  portion  of  the  load/end-shortening  curve  in  the  vicinity  of  the 
c(  uapse  point  is  illustrated  and  the  solutions  obtained  using  the  Riks/Wempner  algorithm  and  the 
ri^rihal  flow  algorithm  (m  =  4)  are  compared.  The  normal  flow  algorithm  was  able  to  traverse  this 
portion  of  the  curve  using  a  relatively  small  number  of  steps  and  continue  on  with  the  remainder 
of  the  curve.  In  contrast,  the  Riks/Wempner  algorithm  was  able  to  make  forward  progress  only 
by  taking  significantly  smaller  steps.  Furthermore,  the  Riks/Wempner  algorithm  was  unable  to 
continue  tracing  the  curve  past  a  certain  point;  at  this  point  the  algorithm  began  chattering  back 
and  forth  over  already  computed  portions  of  the  curve,  perhaps  due  to  the  step  size  control  rules. 
A  Crisfield  solution  obtained  using  equation  (6)  (m  =  4)  and  the  normal  flow  solution  on  the 
same  portion  of  the  curve  are  compared  in  Figure  4.  The  Crisfield  algorithm  is  able  to  compute 
more  of  the  curve  then  the  Riks/Wempner  algorithm,  but  it  also  failed  to  make  forward  progress 
after  reaching  a  certain  point  on  the  curve.  As  with  the  Riks/Wempner  algorithm,  it  is  forced  to 
take  significantly  smaller  (and  therefore  more)  steps,  as  compared  to  the  normal  flow  algorithm,  in 
order  to  obtain  converged  solutions.  Equation  (5)  was  also  used  to  obtain  solutions  the  problem, 
but  this  algorithm  performed  no  better  than  the  equation  (6)  algorithm. 

Using  the  normal  flow  algorithm,  only  one  or  two  iterations  were  required  to  obtain  a  converged 
solution  at  each  step,  whereas  the  Riks/Wempner  and  Crisfield  algorithms  typically  required  from 
three  to  five  iterations  at  each  step.  Depending  on  the  value  of  m,  the  Riks/Wempner  and  Crisfield 
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Figure  4.  Crisfield  trace  of  equilibrium  curve. 

algorithms  each  required  from  250  to  400  matrix  formations/factorizations  in  tracing  the  load/end¬ 
shortening  curve  up  to  the  points  where  the  algorithms  failed;  the  normal  flow  algorithm  reached 
the  same  portion  of  the  curve  with  only  78  formations/factorizations. 

The  increased  computational  efficiency  that  can  be  achieved  using  normal  flow  algorithm  can 
be  particularly  advantagv^uus  in  an  optimization  environment,  where  the  computational  eflBciency 
of  a  single  analysis  is  of  utmost  importance.  If  the  equilibrium  path  must  be  traced  with  a  higher 
degree  of  fidelity  (compared  to  that  shown  in  Figures  3  and  4),  the  normal  flow  algorithm  can 
be  forced  to  take  smaller  steps  along  the  path.  This  will,  of. course,  result  in  a  degradation  of 
computational  efficiency. 

Genetic  algorithms  for  local  panel  design;  The  basic  genetic  algorithm  designed  earlier  in 
the  research  is  capable  of  dealing  with  continuous  design  variables  as  well  as  discrete  variables.  Most 
composite  structural  design  problems  deal  with  determining  the  laminate  stacking  sequences,  which 
is  a  discrete  problem.  However,  continuous  variables,  such  as  stiffener  height  and  core  thickness 
for  sandwich  laminates,  are  often  important  elements  of  a  design  problem.  Therefore,  inclusion 
of  the  continuous  design  variables  into  the  algorithm  was  implemented  by  using  a  state-of-the-art 
approach;  the  continuous  variables  are  handled  as  real  numbers  rather  than  discretizing  them  into 
binary  (or  any  other)  alphabet  strings.  In  addition,  special  operators  for  crossover  and  mutation 
of  the  continuous  numbers  were  used  to  improve  the  efficiency  of  handling  continuous  variables. 

Despite  the  improved  treatment  of  the  continuous  variables,  recent  studies  suggested  that 
the  efficiency  of  the  genetic  algorithm  might  be  increased  for  certain  problems  that  involve  mixed 
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discrete/continuous  variables.  This  is  especially  true  when  there  are  a  large  number  of  near  optimal 
designs  in  the  vicinity  of  the  global  optimal  design.  This  was  the  case  for  a  sandwich  composite 
laminate  design  that  was  performed  for  Sikorsky  Aircraft.  The  problem  was  to  determine  the 
stacking  sequence  of  the  face-sheets  as  well  as  the  thickness  of  the  honeycomb  core  of  the  laminate. 
It  was  observed  that  a  certain  face-sheet  laminate  design  appeared  repeatedly  with  different  values 
of  the  continuous  core  thickness.  In  order  to  improve  the  efficiency  of  the  algorithm,  we  decided  to 
introduce  local  improvements  into  the  search  in  the  form  of  local  approximations  to  the  response  of 
the  panels.  The  approach  used  is  to  monitor  and  ‘’•ore  the  analysis  information  during  the  search 
into  a  binary  tree  structure,  which  serves  as  memo^’y  for  the  optimization. 

For  a  purely  discrete  variable  formulation,  a  binary  tree  memory  structure  is  easy  to  implement, 
and  has  been  used  earlier  by  the  present  investigators  to  improve  the  efficiency  of  search  as  much 
as  30%.  The  improvement  is  due  to  the  fact  that  once  the  analysis  is  performed  for  a  specified 
discrete  combination  of  design  variables,  the  next  time  the  same  combination  appears  the  analysis 
information  can  be  obtained  from  the  memory  instead  of  performing  the  analysis.  In  the  case  of 
mixed  problems,  storing  the  analysis  information  for  every  value  of  the  continuous  design  variables 
would  make  little  sense,  because  even  slight  perturbations  in  the  continuous  values  would  make  a 
design  appear  like  one  that  has  never  been  evaluated.  However,  by  monitoring  the  distribution  and 
proximity  of  the  continuous  design  variables  to  previously  evaluated  -values,  we  should  be  able  to 
construct  local  low  order  polynomial  approximations.  We  are  currently  in  the  process  of  designing 
a  procedure  to  implement  memory  for  the  mixed  discrete/continuous  design  variable  problems. 

Accomplishments/new  findings. 

Genetic  algorithms  for  local  panel  design;  The  basic  genetic  algorithm  data  structures 
were  extended  to  use  multiple  chromosomes  (for  structures  with  several  different  laminates  such  as 
wing  boxes)  and  real  variables  (such  as  geometrical  dimensions) .  The  approach  is  implemented  in 
the  new  object-oriented  Fortran  90  code  for  the  design  of  laminated  composite  panels.  The  multiple 
chromosome  approach  yields  efficient  and  elegant  code,  and  implementing  crossover  and  mutation 
directly  with  real  numbers  is  more  effective  than  ‘i/orking  with  discrete  binary  approximations  of 
real  numbers.  Fortran  90  modules  and  GA  ope-ators  continue  to  be  refined  and  added  to  the 
package,  and  new  capabilities  such  as  memory  (storing  old  designs  in  a  balanced  binary  tree  to 
avoid  reanalysis)  and  local  improvement  (using  local  approximations  to  replace  an  individual  by 
the  best  nearby  individual)  are  being  added. 

Nonlinear  local  panel  analysis:  The  relative  efficiencies  of  the  Riks/Wempner,  Crisfield, 
and  normal  flow  solution  algorithms  for  tracking  nonlinear  equilibrium  paths  of  structural  systems 
are  compared.  It  is  found  that  the  normal  flow  algorithm  may  be  both  more  computationally 
efficient  and  more  robust  compared  to  the  other  two  algorithms  when  tracing  the  path  through 
severe  nonlinearities  such  as  those  associated  with  structural  collapse.  This  is  demonstrated  by 
comparing  the  relative  behaviors  of  each  algorithm  in  the  vicinity  of  a  severe  nonlinearity.  Results 
for  computing  the  collapse  load  of  a  blade  stiffened  panel  demonstrated  substantial  computational 
time  savings  during  the  local  design  of  those  panels. 

Parallel  genetic  algorithms;  The  parallel  GA  code  implemented  migration  between  sub¬ 
populations  (evolving  independently  on  different  processors),  and  the  parallel,  dynamically  adapt¬ 
ing  algorithm  resulted  in  both  better  reliability  and  reduced  normalized  cost  compared  to  a  static 
serial  algorithm  with  migration.  The  parallel  implementation  used  dynamic  load  balancing,  fully 
distributed  control,  and  a  sophisticated  termination  detection  algorithm.  The  surprising  finding 
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was  that  nondeterminism  in  the  parallel  task  management  significantly  enhanced  the  performance 
of  the  (already  nondeterministic)  GA  and  migration.  One  interpretation  of  this  finding  is  that 
random  migration  is  superior  to  a  fixed  migration  pattern.  The  parallel  code  also  scaled  very  well, 
showing  no  significant  communication  penalty  with  64  processors.  Current  work  involves  improv¬ 
ing  the  parallel  GA  code,  exploring  the  effect  of  different  migration  patterns  (corresponding  to 
different  parallel  connection  topologies),  and  blending  local  panel  designs  to  meet  manufacturing 
constraints. 
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Objectives. 


Design  and  optimization  of  large  wing  and  fuselage  structures  presents  a  great  challenge  to 
structural  designers.  Complex  and  highly  local  details  in  these  structures  preclude  the  use  of  a  single 
level  design  model  for  the  entire  structure.  We  have  successfully  introduced  a  two-level  global-local 
design  methodology  to  provide  designers  with  the  capability  of  accounting  for  local  details  (panel 
geometric  variables,  stacking  sequence)  associated  with  local  configurations  (laminated  composite 
sandwich  panels).  At  the  local  level  of  the  two-level  optimization  approach  we  employ  genetic 
algorithms  to  tailor  the  panel  designs.  Response  surface  approximations  to  optimized  panel  failure 
loads  are  then  used  for  the  upper  level  wing  or  fuselage  optimization.  Much  of  the  effort  during  the 
present  reporting  period  was  invested  in  the  coordination  between  codes  that  optimize  the  entire 
wing  or  fuselage  structure  and  codes  that  optimize  individual  panel  details,  so  that  optimal  or  near 
optimal  designs  could  be  obtained.  In  particular,  we  introduced  new  schemes  for  blending  of  local 
laminate  stacking  sequences  across  multiple  panels  using  parallel  computing  and  high-level  object- 
oriented  code  development.  We  also  improved  the  efficiency  of  the  genetic  algorithm  for  local  panel 
designs  by  introducing  a  new  local  improvement  strategy  capable  of  handling  continuous  design 
variables  without  discretization.  Also,  development  of  a  novel  kind  of  design  approach  based  on 
cellular  automata  was  initiated. 

Status  of  research. 

Global/local  optimization  and  blending  using  genetic  algorithms  with  migration: 

The  present  research  concerns  the  design  of  composite  material  aircraft  structures,  particularly 
designing  the  local  details  of  stiffened  and  sandwich  panels  in  the  wing  and  fuselage.  These  panels 
are  made  of  composite  laminates,  using  fiber-reinforced  layers.  The  design  of  the  fiber-reinforced 
composite  laminate  requires  the  specification  of  the  stacking  sequence  with  orientation  and  material 
type,  creating  a  discrete  optimization  problem.  It  is  computationally  expensive  to  design  an  entire 
wing  or  fuselage  structure  with  the  panels  optimized  simultaneously.  Instead,  local  panels  are 
commonly  optimized  (failure  load  is  maximized  given  a  fixed  number  of  plies  of  each  orientation) 
for  the  specified  local  loads.  An  improvement  to  this  process,  a  global/local  procedure,  fits  o 
response  surface  approximation  to  the  maximal  failure  load  of  the  local  panel  configuration 
a  function  of  the  local  loads  and  the  number  of  plies  of  each  orientation.  The  response  surface 
approximation  is  then  used  in  a  global  wing  optimization  (with  numbers  of  plies  of  each  orientation 
as  design  variables)  providing  constraints  on  panel  load  carrying  capacity  while  taking  into  account 
the  variation  of  the  local  load  redistribution  as  a  function  of  the  design  changes.  The  global  solution 
(the  load  distribution)  is  then  used  one  final  time  to  construct  locally  optimal  panels  that  match 
the  overall  wing  design. 

The  global/local  optimization  strategy  briefly  described  above  yields  designs  that  are  locally 
optimal.  However,  process  often  results  in  designs  that  have  incompatibilities  between  the  stacking 
sequence  of  the  adjacent  panels  that  form  the  entire  structure.  Such  incompatibilities  often  result  in 
structures  that  are  difficult  to  manufacture,  or  in  some  cases  impossible  to  manufacture.  Although 
the  differences  are  driven  by  local  loads  and  constraints  that  are  applied  across  the  structure  at 
different  panels  differently,  it  may  be  possible  to  improve  the  manufacturability  by  directly,  or 
indirectly,  coercing  the  adjacent  laminates  to  have  similar  stacking  sequences.  In  particular,  when 
genetic  algorithms  are  used  to  optimize  the  local  panel  stacking  sequences,  often  multiple  near 
optimal  solutions  with  minor  differences  in  the  laminate  stacking  sequences  are  obtained  for  each 
panel.  Therefore,  it  may  be  possible  to  accept  small  losses  from  optimality  of  the  adjacent  panels 
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but  achieve  substantial  gains  in  the  continuity  of  the  laminate  stacking  sequences  between  adjacent 
panels.  In  implementing  the  genetic  algorithm  optimization,  by  allowing  local  designs  to  evolve  in 
parallel,  sending  migrants  to  adjacent  panels,  the  objective  functions  may  be  intelligently  modified 
to  favor  a  globally  blended  design. 

During  the  present  effort  the  blending  of  the  laminate  stacking  sequence  is  accomplished 
using  the  string  edit  distance  between  individuals  and  the  set  of  emigrants  from  adjacent  panels. 
That  is,  instead  of  optimizing  local  panels  in  isolation,  each  individual  is  designed  by  a  different 
piocessor  in  a  parallel  enviroiunent.  The  populations  representing  each  panel  evolve  in  parallel 
and  periodically  send  migrant  individuals  to  adjacent  populations.  A  metric  is  applied,  in  this 
case  the  string  edit  distance,  to  evaluate  the  similarities  in  the  current  population  with  members 
stored  as  migrants.  The  computed  distance  value  is  used  to  reward  (increase  the  fitness  of)  those 
individuals  in  a  population  that  are  evolving  closely  to  those  found  in  the  migrant  populations  from 
adjacent  panels.  This  rewarding  process  is  a  small  modification  to  the  objective  function  fitness 
value  and  serves  to  blend  the  panel  designs  of  the  population  balancing  the  optimization  of  weight 
with  stacking  sequence  compatibility. 

Early  tests  this  year  showed  that  applying  the  blending  metric  to  multiple  migrants  from  a 
single  neighbor  created  a  theoretical  deadlock  situation.  Although  it  is  necessary  to  use  only  a 
single  migrant  from  neighboring  populations  in  the  blending  metric  to  avoid  possible  deadlock, 
this  reduces  genetic  diversity  and  often  results  in  convergence  to  an  inferior  local  minimum.  The 
first  barrage  of  tests  used  an  unmodified  version  of  the  edit  distance  metric.  No  algorithmic 
modifications  were  made  to  match  common  heuristics  used  in  evaluating  the  quality  of  a  blended 
design.  This  would  include  adding  machinery  that  could  evaluate  an  internal  mismatch  in  ply 
angles  as  a  worse  design  than  an  exterior  mismatch.  Similarly  it  could  include  the  matching  of 
an  empty  ply  as  more  desirable  than  a  physical  ply  angle  disruption.  Promising  results  for  the 
method  were  gathered  on  both  two-panel  and  four-panel  ring  topology  test  examples. 

At  this  point  well  blended  designs  became  much  heavier  as  required  plies  were  added  to  satisfy 
the  blending  requirements.  To  combat  this  trend,  heuristic  machinery  was  added  to  the  blending 
algorithm  to  allow  empty  ply  encodings  to  perfectly  match  any  ply  angle  encoding.  A  more 
^'uiplex  and  irregular  grid  example  with  eighteen  panels  was  adopted  for  more  rigorous  testing  qf 
the  methodology  and  is  shown  in  Figure  1.  Ignoring  the  topology  at  first,  tests  were  run  on  the 
eighteen  panels  in  a  ring  formation  with  again  promising  results  and  a  convergence  to  perfectly 
blended  panels  using  the  naive  edit  distance  metric.  Using  the  modified  metric,  allowing  empty 
plies  to  match  anything,  promising  results  were  gained  even  in  a  ring  topology,  but  by  ignoring  the 
true  topology  in  this  way  a  perfectly  blended  design  could  not  be  derived. 

Current  tests  are  running  to  evaluate  a  communication  structure  mapped  directly  from  the 
physical  topology  of  the  panel  layout.  Initial  runs  show  a  much  larger  iteration  count  towards 
convergence  but  should  be  able  to  converge  to  a  perfectly  blended  design  with  respect  to  the 
actual  grid  panel  connections. 

Geometrically  nonlinezu:  analysis  of  local  panels:  During  local  panel  optimization, 
sometimes  it  is  necessary  to  compute  the  nonlinear  relationship  between  the  applied  loads  and 
the  resulting  structural  deformations.  By  tracking  this  nonlinear  equilibrium  path,  it  is  possi¬ 
ble  to  identify  critical  buckling  and  collapse  phenomena.  In  structural  mechanics,  the  nonlinear 
equilibrium  equations  are  often  solved  using  the  Newton  Raphson  method  in  conjunction  with  an 
incremental/iterative  solution  method.  In  order  to  trace  the  equilibrium  paths  through  limit  points 
in  the  response,  the  “arc  length”  methods  of  Riks/Wempner  or  Crisfield  are  commonly  employed. 
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Figure  1.  18  panel  structure  for  testing  parallel  GA  with  blending. 


A  closely  related  method,  which  is  based  on  a  “normal  flow”  algorithm  (from  the  HOMPACK  soft¬ 
ware  package  by  Watson,  Billups,  and  Morgan),  is  less  popular  despite  the  fact  that  it  may  have 
advantages  in  terms  of  computational  efficiency  and  robustness.  During  the  previous  reporting  pe¬ 
riod,  an  implementation  of  the  normal  flow  algorithm  to  nonlinear  local  composite  panel  analyses 
was  investigated,  and  the  potential  advantages  of  the  method  for  the  solution  of  structural  me¬ 
chanics  problems  was  demonstrated.  The  normal  flow  algorithm  for  tracing  nonlinear  equilibrium 
paths,  developed  under  tuic  project,  is  now  in  production  use  for  composite  panel  analysis  in  other 
projects. 

Genetic  algorithms  for  local  panel  design:  The  basic  genetic  algorithm  designed  earlier  in 
the  research  is  capable  of  dealing  with  continuous  design  variables  as  well  as  discrete  variables.  Most 
composite  structural  design  problems  deal  with  determining  thelaminate  stacking  sequences,  which 
is  a  discrete  problem.  However,  continuous  variables,  such  as  stiffener  height  and  core  thickness 
for  sandwich  laminates,  are  often  important  elements  of  a  design  problem.  Therefore,  inclusion 
of  the  continuous  design  variables  into  the  algorithm  was  implemented  by  using  a  state-of-the-art 
approach;  the  continuous  variables  are  handled  as  real  numbers  rather  than  discretizing  them  into 
binary  (or  any  other)  alphabet  strings.  In  addition,  special  operators  for  crossover  and  mutation 
of  the  continuous  numbers  were  used  to  improve  the  efficiency  of  handling  continuous  variables. 

Despite  the  improved  treatment  of  the  continuous  variables,  recent  studies  suggested  that 
the  efficiency  of  the  genetic  algorithm  might  be  increased  for  certain  problems  that  involve  mixed 
discrete/continuous  variables.  This  is  especially  true  when  there  are  a  large  number  of  near  optimal 
designs  in  the  vicinity  of  the  global  optimal  design.  This  was  the  case  for  a  sandwich  composite 
laminate  design  that  was  performed  for  Sikorsky  Aircraft.  The  problem  was  to  determine  the 
stacking  sequence  of  the  face-sheets  as  well  as  the  thickness  of  the  honeycomb  core  of  the  laminate. 
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It  was  observed  that  a  certain  face-sheet  laminate  design  appeared  repeatedly  with  different  values 
of  the  continuous  core  thickness.  In  order  to  improve  the  efficiency  of  the  algorithm,  we  decided  to 
introduce  local  improvements  into  the  search  in  the  form  of  local  approximations  to  the  response  of 
the  panels.  The  approach  used  is  to  monitor  and  store  the  analysis  information  during  the  search 
into  a  binary  tree  structure,  which  serves  as  memory  for  the  optimization. 

For  a  purely  discrete  variable  formulation,  a  binary  tree  memory  structure  is  easy  to  imple¬ 
ment,  and  has  been  used  earlier  by  the  present  investigators  to  improve  the  efficiency  of  search 
as  much  as  30%.  The  improvement  is  due  to  the  fact  that  once  the  analysis  is  performed  for  a 
specified  discrete  combination  of  design  variables.  <^he  next  time  the  same  combination  appears  the 
analysis  information  can  be  obtained  from  the  memory  instead  of  performing  the  analysis.  In  the 
case  of  mixed  problems  with  discrete  and  continuous  variables,  storing  the  analysis  information 
for  every  value  of  the  continuous  design  variables  would  make  little  sense,  because  even  slight 
perturbations  in  the  continuous  values  would  make  a  design  appear  like  one  that  has  never  been 
evaluated.  However,  by  monitoring  the  distribution  and  proximity  of  the  continuous  design  vari¬ 
ables  to  previously  evaluated  values,  accurate  piecewise  polynomial  approximations  to  the  analysis 
dat.a  can  be  constructed.  An  algorithm  to  effect  memory  for  mixed  discrete/continuous  design 
variable  problems  follows. 

Let  n  be  a  discrete  vector  variable,  x  a  scalar  real  variable,  /(v,  x)  the  fitness  of  the  individual 
defined  by  {v,  x),  and  d  a  real  nmnber.  Each  node  in  the  binary  tree  memory  structure  records  the 
tuple  (t;,a;, /(ti,x),d). 

Consider  a  candidate  individual  (v,  x)  for  insertion  in  the  tree. 

If  V  is  not  found  in  the  tree,  then 
evaluate  f{v,  x); 

add  a  node  corresponding  to  {v,  x,  f{v,x),0)] 
return  with  value  f{v,  x)  to  the  GA; 
else 

let  {(n,  Xi,  fi,  be  all  the  node  data  matching  v; 

compute  a  spline  S{-)  of  order  <  4  interpolatmg  the  data  (xj,  fi),  i  —  1,  . . . ,  n; 
evaluate  the  spline  •S'(x); 
if  max  {di  —  |x  —  Xi|}  >0,  then 

return  with  value  S{x)  to  the  GA; 
else 

evaluate /(?;,  x); 
if|/(n,x)  —  5(x)|  >  e,  then 

add  a  node  corresponding  to  {v,  x,  f{v,  x),  0); 
return  with  value  f{v,  x)  to  the  GA; 
else 

define  k  and  d  by  d  =  |x  —  x^l  =  min  |x  —  Xi|; 

l<i<n 

change  the  data  at  node  (v,  Xk,  fk,  dk)  to  {v,  Xk,  fk,  d)-, 
add  a  node  corresponding  to  {v,  x,  f{v,  x),  d); 
return  with  value  f{v,  x)  to  the  GA; 
end  if 
end  if 
end  if 
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During  the  last  activity  period,  an  additional  use  of  the  spline  approximation  for  the  continuous 
variables  was  developed.  Rather  then  using  the  spline  passively  just  as  a  memory  (or  surrogate 
analysis)  for  new  design  points,  we  have  started  using  it  actively  in  generating  new  points  that 
will  be  included  into  the  next  generation.  That  is,  knowing  the  spline  curve  in  a  given  generation 
we  were  able  to  determine  the  value  of  the  continuous  variable  that  provided  the  best  value  of  the 
fitness  function  for  a  particular  value  of  the  discrete  variable  associated  with  the  spline  (note  that 
there  is  a  spline  approximation  for  each  discrete  node  in  the  binary  tree  if  that  node  is  visited 
more  than  twice).  When  creating  children  from  the  parent  generation  ’r-iing  cross-over,  if  the 
discrete  part  of  the  child  was  associated  with  a  spline  that  stored  the  best  value  of  the  continuous 
variable,  then  this  value  is  used  to  replace  the  value  of  the  continuous  cross-over  value  for  the 
child  design.  Results  generated  so  far  indicate  substantial  improvement  in  the  performance  of  the 
genetic  algorithm  in  terms  of  the  number  of  function  evaluations  to  reach  an  optimal  design. 

Structural  Design  Using  Cellular  Automata:  During  the  present  performance  period, 
a  new  task  was  initiated.  While  the  traditional  method  of  structural  design  using  finite  element 
based  numerical  analysis  programs  works  well  for  many  problems,  it  does  not  parallelize  efficiently 
on  massively  parallel  processors  (MPPs),  thus  limiting  the  size  and  complexity  of  the  structures 
that  can  be  designed.  The  new  task  demonstrates  the  use  of  a  new  implementation  of  the  cellular 
computational  paradigm  for  simultaneous  analysis  and  design  of,  two-dimensional  structural  do¬ 
mains.  In  the  following,  a  brief  description  of  the  method  and  its  adaptation  to  structural  analysis 
and  design  is  provided.  An  example  is  provided,  and  its  convergence  properties  are  discussed. 

Method  Description:  Cellular  automata  (CA)  were  used  at  least  as  early  as  1946  by  Wiener  [4] 
to  describe  the  operation  of  heart  muscle,  even  though  their  use  was  not  computationally  feasible 
at  the  time.  CA  tiles  a  problem  domain  into  cells  of  equal  size.  Each  cell  has  the  same  set  of 
simple  rules  that  dictate  how  it  behaves  and  interacts  with  its  neighboring  cells.  Each  cell  is  a 
fixed  point  in  a  regular  lattice.  The  state  of  each  cell  is  updated  at  discrete  time  steps,  based  upon 
conditions  in  previous  time  steps.  All  of  the  cells  are  updated  every  time  step,  thus  the  state  of 
the  entire  lattice  is  updated  every  time  step.  The  basic  principle  behind  the  method  is  that  an 
overall  global  behavior  can  be  computed  by  a  group  of  cells  that  only  knov.iocal  conditions  [5].  If 
each  cell  only  needs  to  know  local  conditions,  then  this  minimizes  the  communication  requirements 
and  therefore  the  problem  scales  well  on  a  MPR  A  CA  is  the  archetypical  algorithm  for  the  SIMD 
parallel  architecture  [3]. 

CA  have  recently  been  used  to  perform  structural  analysis  and  design  [1,  2].  In  this  particular 
implementation,  the  intention  is  to  describe  the  static  equilibrium  of  a  two-dimensional  domain 
under  a  system  of  forces  acting  on  it.  In  this  sense,  time  is  not  being  simulated,  rather  each  step 
of  the  automaton  is  used  to  propagate  (local)  stresses  and  strains  through  out  the  domain  to  allow 
it  to  reach  equilibrium  state  while  simultaneously  determining  the  shape  and/or  dimensions  of  the 
cells  associated  with  this  equilibrium  state.  This  is  continued  until  the  entire  process  converges 
(ideally)  to  a  global  state  where  there  is  no  significant  change  in  the  structure  for  every  subsequent 
iteration,  corresponding  to  a  static  equilibrium  state.  Note  that,  as  will  be  described  later,  analysis 
and  design  are  done  simultaneously  and  locally  by  each  cell. 

Domain  Definition:  Each  cell  of  this  CA  is  an  eight-beam  truss  where  each  beam  starts  at  the 
center  of  the  cell  and  connects  to  its  opposite  member  in  an  adjacent  cell  as  illustrated  in  Figure 
2.  This  t5q)e  of  structure  is  known  as  a  ground  truss.  Those  cells  that  fall  on  the  border  of  the 
rectangular  domain  are  not  partial  cells  requiring  special  rules,  but  are  complete  cells  with  the 
area  of  the  beams  that  fall  outside  the  computational  domain  set  to  zero.  In  addition,  they  are 
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Figure  2.  Example  CA  domain  with  a  single  cell  denoted  by  the  dashed  line. 

connected  to  an  invisible  set  of  surrounding  cells  that  are  turned  off  and  that  also  have  all  of  their 
beam  areas  set  to  zero.  Cells  that  are  turned  off  are  not  part  of  the  computation,  being  used  only 
to  make  the  rules  for  the  border  and  non-border  cells  consistent,  since  the  stress  analysis  rules 
require  the  displacements  of  all  eight  surrounding  cells. 

The  actual  border  of  the  computational  domain  of  the  CA  need  not  be  rectangular.  Any  shape 
can  be  defined  for  the  truss  by  turning  off  any  cells  that  are  not  within  the  computational  domain, 
as  illustrated  in  Figure  2.  A  simple  method  to  define  a  shape  for  the  truss  is  to  define  an  enclosing 
polygon,  and  then  turn  off  every  cell  that  does  not  fall  within  the  polygon.  The  “edge  crossings’" 
algorithm  to  determine  those  points  within  the  polygon  can  be  used;  it  is  simple  and  parallelism 
well.  A  more  sophisticated  method  could  be  used  to  allow  for  holes,  circular  regions,  or  other, 
nonstraight  boundaries.  In  addition,  a  better  resolution  can  be  obtained  by  decreasing  the  cell  size 
in  the  domain,  thereby  increasing  the  number  of  cells  that  form  the  shape.  The  amount  by  which 
the  original  cells  have  been  subdivided  to  increase  the  resolution  is  known  as  the  cell  density  factor 
(CDF). 

CA  Rules:  There  are  two  sets  of  rules  employed  to  compute  the  response  of  the  truss  to  applied 
loads  or  displacements,  and  the  change  in  the  cross  sectional  areas  of  the  truss  members  to  fully 
utilize  the  members.  The  first  set  of  rules  is  (normally)  executed  at  every  iteration  to  determine 
the  displacements  of  a  cell  eissociated  with  deformations.  The  cell  attempts  to  reach  equilibrium 
with  the  surrounding  cells  by  displacing  itself  to  minimize  the  potential  energy. 

Within  a  cell,  each  truss  member  (indexed  relative  to  the  cell  by  A:  =  1, . . . ,  8)  has  an  elastic 
modulus  E,  length  Lk,  a  cross-sectional  area  Ak,  and  an  orientation  angle  9k  from  the  cell  center. 
Denote  the  displacement  of  the  fe-th  truss  member’s  near  end  from  the  original  cell  center  by  u,  v, 
and  the  displacement  of  the  neighboring  cell’s  center  by  Uk,Vk.  These  neighboring  displacements 
are  assumed  to  be  specified  when  the  CA  calculates  the  displacements  of  a  cell.  The  extension  of 
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a  member  between  two  cells  Afc,  strain  ejt,  and  force  Fk  within  each  member  are  calculated  from 
these  properties  and  displacements  by 


Afc  =  (u  -  Uc)  cosOk  +  (u  -  Vc)  sin0fc. 

(1) 

=  ^k/Lki 

(2) 

Fk  =  EAkSk- 

(3) 

Taking  into  account  the  applied  external  forces  Fj;,Fy,  the  total  (internal  strain  plus  external) 
potential  energy  V  for  a  cell  is  given  by 

8 

V  =  Y^  EAkLksl  -F^{u-  Uc)  -Fy{v-  v^). 
k=l 

Setting  the  partial  derivatives  of  the  potential  energy  with  respect  to  the  cell  displacements  to  zero 
gives  the  equilibrium  equations  that  can  be  solved  for  the  cell  displacements.  In  general,  this  is  a 
system  of  two  equations  with  two  unknowns.  If  there  is  an  (externally)  applied  displacement  along 
one  or  both  of  the  axes,  then  the  equations  will  be  solved  for  forces,  Fx  and/or  Fy,  instead  of  the 
displacements. 

Designing  the  structure  requires  resizing  the  beams  in  the  cells.  If  displacements  have  already 
been  calculated  then  some  scheme  for  changing  the  cross  sectional  areas  Ak  is  required.  In  terms 
of  allowable  stress  cTaiiow,  which  is  specified  as  the  maximum  stress  that  any  given  beam  should 
endure,  one  scheme  for  computing  a  new  cross  sectional  area  ,  based  upon  the  previous  cross 
sectional  area  A'j^'^,  is 

/inew  _  E\ek\  A  old 

-  ~  • 

C^allow 

If  the  displacement  calculation  and  sizing  are  done  sequentially,  the  sizing  period  (how  often  sizing 
done)  depends  on  many  factors:  the  number  of  cells  in  the  domain,  the  locations  and  relative 
r  lo^cements  of  the  applied  forces  and  displacements,  the  iteration  method  (Jacobi  vs.  Gauss-Seidel), 
and  for  Gauss-Seidel  implemented  in  parallel,  the  number  of  processors  used. 

Example:  Consider  the  problem  of  a  simple  bridge  truss.  The  first  image  in  Figure  3  shows  a 
CA  with  six  cells.  The  bottom  two  corner  cells  have  an  applied  displacement  of  (0,0)  so  they  are 
fixed  in  place.  The  bottom  middle  cell  has  an  applied  force  of  lOOkN  downward.  The  width  of 
the  bridge  is  50  meters  and  the  height  is  25  meters.  The  bars  are  composed  of  medium  steel  with 
E  =  200GPa  and  cTaiiow  =  250MPa.  Each  beam  has  an  initial  area  of  0.0175m^. 

Running  the  CA  on  the  bridge  problem  using  the  Gauss-Seidel  iteration  method  for  displace¬ 
ments  and  applying  the  sizing  rules  every  sixth  iteration  until  it  converges  at  iteration  253,  the 
result  shown  in  the  second  image  of  Figure  3  is  obtained.  Since  the  bridge  is  50m  across  and  the 
steel  beams  are  no  more  than  a  few  cm  thick,  the  areas  in  this  view  are  exaggerated  by  a  factor  of 
3000  to  show  the  differences  in  the  beam  sizes. 

The  bridge  in  Figure  3  is  only  composed  of  6  cells,  and  the  solution  could  easily  have  been 
computed  by  hand.  By  increasing  the  cell  density,  the  complexity  of  the  problem  rises.  Figure 
4  shows  the  final  configuration  of  a  problem  of  the  exact  same  dimensions,  where  each  cell  is  40 
times  smaller  than  previously.  Each  horizontal  and  vertical  beam  is  0.625m  long,  rather  than  25m. 
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Figure  3.  Simple  bridge  truss,  before  and  after  running  CA. 


Figure  4.  Bridge  truss  with  3321  cells  after  convergence. 

Convergence  Analysis:  To  analyze  the  efficacy  of  the  iteration  method  used,  it  is  useful  to 
transform  the  CA  into  an  equivalent  system  of  linear  equations.  Recalling  that  each  cell  is  com¬ 
puting  its  position  n,  v  based  upon  the  position  of  the  surrounding  cells.  If  each  cell  were  assigned 
unique  variables  for  its  position,  such  that  cell  1  has  ui  and  vi ,  cell  2  has  U2  and  V2 ,  and  so  forth, 
then  the  equations  for  each  cell  can  be  expressed  in  terms  of  the  variables  for  the  surrounding 
cells.  For  a  CA  structure  composed  of  6  cells,  this  will  form  a  linear  system  of  12  equations  and 
12  unknowns.  This  standard  system  of  linear  equations,  Ax  =  b,  can  be  solved  by  the  Jacobi  and 
Gauss-Seidel  fixed-point  iteration  methods  or  block  versions  thereof,  which  are  the  exact  mathe¬ 
matical  formulations  of  the  local  cell  calculations.  For  the  Jacobi,  A  is  split  into  its  strictly  2x2 
block  lower  triangular  (L),  2  x  2  block  diagonal  (D),  and  strictly  2x2  block  upper  triangular  (U) 
parts, 

A=L+D+U 

The  system  is  then  rewritten  as  a  fixed  point  iteration  where  the  next  iterate  is  computed 

from  the  previous  iterate  via 


x(n+l)  ^  ^  c, 


CDF. 

n 

Jacobi 

Gauss-Seidel 

1 

8 

0.794104 

0.611503 

2 

26 

0.949748 

0.8982?? 

3 

52 

0.979368 

0.959006 

4 

86 

0.989496 

0.978893 

5 

128 

0.993801 

0.987605 

6 

178 

0.995959 

0.991925 

7 

236 

0.997181 

0.994365 

8 

302 

0.99793.^ 

0.995865 

9 

376 

0.998426 

0.996853 

10 

458 

0.998765 

0.997532 

Table  1.  Spectral  radius  of  the  bridge  truss  for  various  CDFs. 

where 

B  =  -D-^(L  +  U),  C  =  D-^b. 

Note  that  Ax  =  b  if  and  only  if  x  =  Bx  +  C,  assuming  exists.  For  Gauss-Seidel  the 
iteration  =  Bx^”^  +  C  has 


B  =  -(D  +  L)-'U,  C  =  (D  +  L)-ib, 


assuming  (D  +  L)  ^  exists. 

The  fixed-point  iteration  converges  for  any  starting  point  x^^^  if  and  only  if  all  of  the  eigenvalues 
of  B  are  less  than  one  in  absolute  value.  The  maximum  absolute  value  of  the  eigenvalues  for  a 
matrix  is  called  the  spectral  radius.  The  spectral  radius  for  the  bridge  structure  at  various  cell 
density  factors  (CDF)s  is  shown  in  Table  1. 

This  table  shows  that  the  Jacobi  or  Gauss-Seidel  CA  iteration  for  analysis  does  converge,  but 
extremely  slowly.  For  larger  CDFs,  the  improven  ^iit  of  Gauss-Seidel  over  Jacobi  is  marginal.  Even 
with  massive  parallelism,  any  competitive  adva..itage  of  CA  (over  solving  the  linear  system  with 
standard  iterative  numerical  methods)  must  come  by  combining  analysis  with  sizing. 

[1]  Z.  Giirdal  and  B.  Tatting,  “Cellular  automata  for  design  of  truss  structures  with  linear  and 

nonlinear  response” ,  AIAA-2000-1580,  AIAA  Structures,  Structural  Dynamics,  and  Materials 
Conference,  April,  2000,  Atlanta,  GA. 

[2]  Z.  Giirdal  and  B.  Tatting,  “Cellular  automata  for  design  of  two-dimensional  continuum  struc¬ 

tures”,  AIAA- 2000-4832,  8th  AIAA/NASA/ISSMO  Symposium  on  Multidisciplinary  Analysis 
and  Optimization,  September  2000,  Long  Beach  CA. 

[3]  T.  Toffoli  and  N.  Margolus,  “Cellular  Automata  Machines”,  MIT  Press,  Cambridge,  Mas¬ 

sachusetts,  1991. 

[4]  N.  Wiener  and  A.  Rosenblunth,  “The  mathematical  formulation  of  the  problem  of  conduction 

of  impulses  in  a  network  of  connected  excitable  elements,  specifically  in  cardiac  muscle” ,  Arch. 
Inst.  Cardiol.  Mexico,  vol.  16,  pp.  205-265,  1946. 

[5]  S.  Wolfram,  “Cellular  Automata  and  Complexity;  Collected  Papers”,  Addison- Wesley  Publish¬ 

ing  Company,  Reading  Massachusetts,  1994. 
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Accomplishments/new  findings. 

Genetic  algorithms  for  local  panel  design:  The  basic  genetic  algorithm  data  structures 
were  extended  to  use  multiple  chromosomes  (for  structures  with  several  different  laminates  such  as 
wing  boxes)  and  real  variables  (such  as  geometrical  dimensions).  The  approach  is  implemented  in 
the  new  object-oriented  Fortran  90  code  for  the  design  of  laminated  composite  panels.  The  multiple 
chromosome  approach  yields  efficient  and  elegant  code,  and  implementing  crossover  and  mutation 
directly  with  real  nmnbers  is  more  effective  than  working  with  discrete  binary  approximations 
of  real  numbers.  Fortran  90  modules  and  GA  operators  continue  to  be  refined  and  added  to  the 
package,  and  new  capabilities  such  as  discrete  memory  (storing  old  designs  in  a  balanced  binary  tree 
to  avoid  reanalysis)  and  continuous  memory  (extending  the  binary  tree  via  spline  approximations 
in  the  continuous  variable)  were  added.  Local  improvement  (using  local  approximations  to  replace 
an  individual  by  the  best  nearby  individual)  is  also  being  added. 

Laminate  blending  with  Genetic  Algorithms:  In  parallel  with  the  research  carried  out 
at  Virginia  Tech,  a  simple  blending  algorithm  that  is  based  on  an  automated  2-step  optimization 
process  was  incorporated  by  ADOPTECH  Inc.  into  a  commercial  genetic  algorithm  software  code 
called  DARWIN.  The  approach  is  based  on  the  observation  that  if  the  individual  panels  of  an 
assembly  of  panels  are  optimally  designed  by  completely  ignoring  the  need  for  blending,  the  total 
weight  of  the  panel  system  will  be  the  lowest  possible  weight.  However,  the  stacking  sequence 
of  the  individually  designed  panels  will  rarely  have  any  continuity  across  the  panels,  potentially 
depicting  the  least  amount  of  blending  characteristics.  Alternatively,  a  single  sub-laminate  could 
be  defined  for  all  panels  without  any  prior  knowledge  about  the  characteristics  of  the  individually 
optimized  designs.  In  this  case  the  panels  will  have  the  best  blending  characteristics,  and  the 
weigh  will  constitute  an  upper  bound  weight  since  the  most  critical  panel  will  dictate  the  laminate 
stacking  sequence  of  the  entire  panel  system.  Between  these  two  extreme  cases,  other  solutions 
exist  that  produce  a  total  weight  closer  to  that  of  the  lower  bound  weight  and  display  good  blending 
characteristics.  The  approach  used  by  ADOPTECH  is  a  systematic  procedure  that  searches  for  a 
blended  system  of  panels  by  starting  from  the  lower  bound  design. 

The  methodology  is  tested  using  a  3  x  3  array  of  sandwich  panels,  and  a’',  18-panel  group  (Figure 
1)  arranged  in  a  horseshoe  pattern  provided  to  ADOPTECH  by  Sikorsky  Aircraft.  Results  for  eack 
design  problem  are  generated  using  ADOPTECH’s  Java-based  composite  laminate  design  software 
called  OLGA,  which  utilizes  the  DARWIN  optimization  engine.  Those  results  are  presented  in  a 
paper  jointly  authored  by  ADOPTECH,  Sikorsky,  and  Virginia  Tech  personnel. 

Nonlinear  local  peuiel  analysis:  Earlier  work  under  this  grant  compared  the  relative  effi¬ 
ciencies  of  the  Riks/Wempner,  Crisfield,  and  normal  flow  solution  algorithms  for  tracking  nonlinear 
equilibrium  paths  of  structural  systems.  It  was  found  that  the  normal  flow  algorithm  was  both 
more  computationally  efficient  and  more  robust  than  the  other  two  algorithms  when  tracing  the 
path  through  severe  nonlinearities  such  as  those  associated  with  structural  collapse  of  a  blade  stiff¬ 
ened  panel.  Consequently  the  normal  flow  algorithm  is  being  routinely  used  in  other  projects  for 
tracking  equilibrium  curves. 

Parallel  genetic  algorithms:  The  parallel  GA  code  implemented  migration  between  sub¬ 
populations  (evolving  independently  on  different  processors),  and  the  parallel,  dynamically  adapt¬ 
ing  algorithm  resulted  in  both  better  reliability  and  reduced  normalized  cost  compared  to  a  static 
serial  algorithm  with  migration.  The  parallel  implementation  used  dynamic  load  balancing,  fully 
distributed  control,  and  a  sophisticated  termination  detection  algorithm.  The  surprising  finding 
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was  that  nondeterminism  in  the  parallel  task  management  significantly  enhanced  the  performance 
of  the  (already  nondeterministic)  GA  and  migration.  One  interpretation  of  this  finding  is  that 
random  migration  is  superior  to  a  fixed  migration  pattern.  The  parallel  code  also  scaled  very  well, 
showing  no  significant  communication  penalty  with  64  processors.  Current  work  involves  improv¬ 
ing  the  parallel  GA  code,  exploring  the  effect  of  different  migration  patterns  (corresponding  to 
different  parallel  connection  topologies),  and  blending  local  panel  designs  to  meet  manufacturing 
constraints. 

Mathematical  theory  for  global/local  decomposition:  A  rigorous  decomposition  theory 
has  been  developed  for  a  class  of  optimization  problems  under  mild  smoothness  assumptions. 
Necessary  and  suflScient  conditions  under  which  the  decomposition  reformulation  is  equivalent  to 
the  original  problem  have  been  derived.  The  reformulation  can  potentially  remove  a  large  number 
of  local  optima  from  the  original  problem,  making  the  decomposition  computationally  attractive. 
The  next  task  is  to  extend  these  results  to  mixed  continuous/discrete  problems,  so  as  to  include 
inherently  discrete  panel  design. 
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Haftka,  “VizCraft:  a  problem  solving  environment  for  configuration  design  of  a  high  speed 
civil  transport”,  Comput.  Sci.  Engrg.,  3  (2001)  56-66. 

M.  Sosonkina,  J.  T.  Melson,  Y.  Saad,  and  L.  T.  Watson,  “Preconditioning  strategies  for  linear 
systems  arising  in  tire  design”,  Numer.  Linear  Algebra  Appl.,  7  (2000)  743-757. 

L.  T.  Watson  and  C.  A.  Baker,  “A  fully-distributed  parallel  global  search  algorithm”,  Engrg. 
Comput.,  18  (2001)  155-169. 
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Refereed  conference  papers  published  during  the  grant  period  are: 

A.  Goel,  C.  Phanouriou,  F.  A.  Kamke,  C.  J.  Ribbens,  C.  A.  Shaffer,  and  L.  T.  Watson,  “WBCSim: 
A  prototype  PSE  for  wood-based  composites  simulations” ,  in  Enabling  Technologies  for  Com¬ 
putational  Science;  Frameworks,  Middleware  and  Environments,  E.  N.  Houstis,  J.  R.  Rice,  E. 
Gallopoulos,  and  R.  Bramley  (eds.),  Kluwer,  Norwell,  MA,  2000,  187-196. 

C.  A.  Shaffer,  N.  Ramakrishnan,  L.  T.  Watson,  R.  Dymond,  and  V.  Lohani,  “The  role  of  prob¬ 

lem  solving  environments  in  watershed  assessment” ,  in  Proc.  Integrated  Decision-Making  for 
Watershed  Management  Symposium:  Processes  and  Tools,  Chevy  Chase,  MD,  2001,  677-690. 

E.  J.  Rubin,  R.  Dietz,  J.  Chanat,  C.  Speir,  R.  Dymond,  V.  Lohani,  D.  Kibler,  D.  Bosch,  C.  A. 
Shaffer,  N.  Ramakrishnan,  and  L.  T.  Watson,  “From  landscapes  to  waterscapes:  a  PSE  for 
landuse  change  analysis” ,  in  Proc.  Integrated  Decision-Making  for  Watershed  Management 
Symposium:  Processes  and  Tools,  Chevy  Chase,  MD,  2001,  643-659. 

D.  B.  Adams,  L.  T.  Watson,  and  Z.  Giirdal,  “Global/local  optimization  and  blending  using  genetic 

algorithms  with  migration”,  in  Proc.  10th  SIAM  Conf  on  Parallel  Processing  for  Scientific 
Computing,  CD-ROM,  SIAM,  Philadelphia,  PA,  2001,  8  pages. 

L.  T.  Watson  and  C.  A.  Baker,  “A  fully  distributed  parallel  global  search  algorithm” ,  in  Proc.  10th 
SIAM  Conf  on  Parallel  Processing  for  Scientific  Computing,  CD-ROM,  SIAM,  Philadelphia, 
PA,  2001,  10  pages. 

J.  W.  Zwolak,  J.  J.  Tyson,  and  L.  T.  Watson,  “Estimating  rate  constants  in  cell  cycle  models” ,  in 
Proc.  High  Performance  Computing  Symposium  2001,  A.  Tentner  (ed.),  Soc.  for  Modeling 
and  Simulation  Internat.,  San  Diego,  CA,  2001,  53-57. 

A.  Verstak,  M.  Vass,  N.  Ramakrishnan,  C.  Shaffer,  L.  T.  Watson,  K.  K.  Bae,  J.  Jiang,  W.  H. 
Tranter,  and  T.  S.  Rappaport,  “Lightweight  data  management  for  compositional  modeling  in 
,  problem  solving  environments”,  in  Proc.  High  Performance  Computing  Symposium  2001,  A. 
Tentner  (ed.),  Soc.  for  Modeling  and  Simulation  Internat.,  San  Diego,  CA,  2001,  148-153. 

Z.  Giirdal  and  L.  T.  Watson,  “Structural  design  using  cellular  automata” ,  in  Proc.  Third  Internat. 
Conf.  on  Intelligent  Processing  and  Manufacturing  of  Materials,  J.  A.  Meech,  S.  M.  Veiga, 
M.  M.  Veiga,  S.  R.  LeClair,  and  J.  F.  Maguire  (eds.),  CD-ROM,  Vancouver,  Canada,  2001,  4 
pages. 

W.  Chauncey,  and  Z.  Gurdal  “Thermal  testing  of  tow-placed  variable  stiffness  panels,”  AIAA- 
2001-1190,  in  Proc.  42nd  AIAA/ASME/ASCE/AHS/ASC  Structures,  Structural  Dynamics, 
and  Materials  Conf,  CD-ROM,  April  16-19,  Seattle,  WA,  2001. 

S.  Missoum  and  Z.  Gurdal  “A  displacement  based  optimization  for  nonlinear  frame  structures,” 
AIAA-2001-1553,  in  Proc.  42nd  AIAA/ASME/ASCE/AHS/ASC  Structures,  Structural  Dy¬ 
namics,  and  Materials  Conf,  CD-ROM,  April  16-19,  Seattle,  WA,  2001. 

G.  Soremekun,  Z.  Gurdal,  C.  Kassapoglou,  and  D.  Toni,  “Stacking  sequence  blending  of  multi¬ 
ple  composite  laminates  using  genetic  algorithms,”  AIAA-2001-1203,  in  Proc.  42nd  AIAA/- 
ASME/ASCE/AHS/ASC  Structures,  Structural  Dynamics,  and  Materials  Conf,  CD-ROM, 
April  16-19,  Seattle,  WA,  2001. 


13 


Journal  articles  accepted  during  the  grant  period  are: 

C.  A.  Baker,  B.  Grossman,  R.  T.  Haftka,  W.  H.  Mason,  and  L.  T.  Watson,  “HSCT  configuration 
design  space  exploration  using  aerodynamic  response  surface  approximations” ,  J.  Aircraft,  to 
appear. 

S.  A.  Ragon,  Z.  Giirdal,  and  L.  T.  Watson,  “A  comparison  of  three  algorithms  for  tracing  nonlinear 
equilibrium  paths  of  structural  systems” ,  Internat.  J.  Solids  Structures,  to  appear. 

L.  T.  Watson,  “Theory  o^  globally  convergent  probability-one  homotopies  for  nonlinear  program¬ 
ming”,  SIAM  J.  Optim.,  to  appear. 

C.  A.  Baker,  L.  T.  Watson,  B.  Grossman,  W.  H.  Mason,  and  R.  T.  Haftka,  “Parallel  global  aircraft 

configuration  design  space  exploration”,  Internat.  J.  Comput.  Res.,  to  appear. 

B.  G.  Zombori,  F.  A.  Kamke,  and  L.  T.  Watson,  “Modeling  the  transient  effects  during  the  hot- 
pressing  of  wood-based  composites — Part  I:  simulation  of  the  mat  formation  process” ,  Wood 
Fiber  Sci.,  to  appear. 

L.  T.  Watson,  “Probability-one  homotopies  in  computational  science” ,  J.  Comput.  Appl.  Math., 

to  appear. 

H.  Kim,  M.  Papila,  W.  H.  Mason,  R.  T.  Haftka,  L.  T.  Watson,  and  B.  Grossman,  “Detection  and 
repair  of  poorly  converged  optimization  runs”,  AIAA  J.,  to  appear. 

S.  C.  Billups  and  L.  T.  Watson,  “A  probability-one  homotopy  algorithm  for  nonsmooth  equations 
and  mixed  complementarity  problems”,  SIAM  J.  Optim.,  to  appear. 

S.  Hosder,  L.  T.  Watson,  B.  Grossman,  W.  H.  Mason,  H.  Kim,  R.  T.  Haftka,  and  S.  Cox,  “Polyno¬ 
mial  response  surface  approximations  for  the  multidisciplinary  design  optimization  of  a  high 
speed  civil  transport”,  Optim.  Engrg.,  to  appear. 

S.  Missoum,  Z.  Giirdal,  and  G.  Weinjiong,  “Optimization  of  nonlinear  trusses  using  a  displacement 
baaed  optimization,”  Structural  and  Multidisciplinary  Optim.,  to  appear. 

S. 'Missoum  and  Z.  Giirdil,  “A  displacement  based  optimization  for  trusses  structures  subjected 
^  to  static  and  dynarric  constraints,”  AIAA  J.,  to  appear. 

Journal  articles  submitted  during  the  grant  period  are: 

M.  Sosonkina,  D.  C.  S.  Allison,  and  L.  T.  Watson,  “Scalability  analysis  of  parallel  GMRES  imple¬ 

mentations”  ,  Parallel  Algorithms  Appl. 

G.  Mateescu,  C.  J.  Ribbens,  and  L.  T.  Watson,  “A  domain  decomposition  preconditioner  for 

Hermite  collocation  problems” ,  Numer.  Methods  Partial  Differential  Equations. 

S.  E.  Cox,  R.  T.  Haftka,  C.  A.  Baker,  B.  Grossman,  W.  H.  Mason,  and  L.  T.  Watson,  “Global 
multidisciplinary  optimization  of  a  high  speed  civil  transport” ,  J.  Global  Optim. 

H.  Kim,  R.  T.  Haftka,  W.  H.  Mason,  L.  T.  Watson,  and  B.  Grossman,  “A  study  of  the  statistical 

modelling  of  errors  from  structural  optimization” ,  Optim.  Engrg. 

E.  J.  Rubin,  R.  Dietz,  J.  Chanat,  C.  Speir,  R.  Dymond,  V.  Lohani,  D.  Kibler,  D.  Bosch,  C.  A. 
Shaffer,  N.  Ramakrishnan,  and  L.  T.  Watson,  “Prom  landscapes  to  waterscapes:  a  PSE  for 
landuse  change  analysis” ,  Engrg.  Computers. 

D.  J.  Slotta,  B.  Tatting,  L.  T.  Watson,  and  Z.  Giirdal,  “Structural  design  using  cellular  automata”, 

Engrg.  Comput. 
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C.  A.  Shaffer,  N.  Ramakrishnan,  L.  T.  Watson,  R.  Dymond,  and  V.  Lohani,  “The  role  of  problem 

solving  environments  in  watershed  assessment” ,  J.  Amer.  Water  Resources  Assoc. 

L.  T.  Watson,  V.  K.  Lohani,  D.  F.  Kibler,  R.  L.  Dymond,  N.  Ramakrishnan,  and  C.  A.  Shaffer, 

“Integrated  computing  environments  for  watershed  management”,  J.  Comput.  Civil  Engrg. 

M.  A.  Gongalves,  E.  A.  Fox,  L.  T.  Watson,  and  N.  A.  Kipp,  “Streams,  structures,  spaces,  scenarios, 

societies  (5S):  a  formal  model  for  digital  libraries”,  ACM  Trans.  Infor.  Systems. 

D.  G.  Kafura,  N.  Ramakrishnan,  C.  J.  Ribbens,  G,  A.  Shaffer,  and  L.  T.  Watson,  “Programming 

environments  for  multidisciplinary  grid  communities” ,  Concurrency:  Pract.  Exper. 

S.  Missoum,  Z.  Giirdal,  and  L.  T.  Watson,  “A  displacement  based  optimization  method  for  geo¬ 
metrically  nonlinear  frame  structures” ,  Structural  Multidisciplinary  Optim. 

G.  Weinjiong,  Z.  Giirdal,  and  S.  Missoum,  “Holonomic  elastoplastic  truss  design  using  displacement 
based  optimization,”  J.  Comp.  Meth. 

.  S.  Ragon,  S.  Chandrasekaran,  D.  K.  Lindner,  Z.  Giirdal,  and  D.  Boroyevich,  “Optimal  design  of 
a  power  distribution  system,”  J.  Aircraft. 

Interactions  /  transitions. 

Technology  transitions  or  transfer: 

PERFORMER 

Layne  T.  Watson,  Virginia  Polytechnic  Institute  &  State  University 

Telephone:  540-231-7540 

CUSTOMER 

General  Motors  Research  and  Development  Center 
Warren,  MI 

Contact:  Alexander  P.  Morgan,  810-986-2157 
RESULT 

Homotopy  algorithms;  mathematical  software 
APPLICATION 

Linkage  mechanism  design;  combustion  chemistry;  robotics;  CAD/CAM 
PERFORMER 

Layne  T.  Watson,  Virginia  Polytechnic  Institute  Sc  State  University 

Telephone:  540-231-7540 

CUSTOMER 

Lucent  Technologies 

Murray  Hill,  NJ 

Contact:  Robert  Melville,  908-582-2420 
RESULT 

Homotopy  algorithms;  mathematical  software 

APPLICATION 

Circuit  design  and  modelling 
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PERFORMER 

Layne  T.  Watson,  Virginia  Polytechnic  Institute  &  State  University 

Telephone:  540-231-7540 

CUSTOMER 

Michelin  Americas 

Greenville,  SC 

Contact:  John  Melson,  864-422-4246 
RESULT 

Adaptive  GMRES  algorithm;  mathematical  software 
APPLICATION 

Iterative  solution  of  large  linear  systems  arising  from  tire  modelling 
PERFORMER 

Zafer  Giirdal,  Virginia  Polytechnic  Institute  &  State  University 

Telephone:  540-231-5905 

CUSTOMER 

Sikorsky  Aircraft 

Stratford,  Connecticut  06497 

Contact:  Christos  Kassapoglou,  203-386-3292 

RESULT 

Multiobjective  genetic  algorithms  for  composite  laminate  design 
APPLICATION 

Design  of  helicopter  frame  structures  for  minimum  weight  and  cost 
PERFORMER 

Zafer  Giirdal,  Virginia  Polytechnic  Institute  &  State  University 
Telephone:  540-231-5905 
CUSTOMER 
ADOPTech  Inc. 

Blacksburg,  VA  24061 

Contact:  Grant  Soremekun,  540-231-7232 

RESULT 

Blending  algorithm  for  multipanel  composite  structural  design 
APPLICATION 

Design  of  helicopter  skin  structures  for  minimum  weight  and  cost 

Inventions  or  patents. 


None. 


Honors  /  awards. 

•  IEEE  Fellow:  Layne  T.  Watson. 

•  Best  Paper  Award,  Sixth  AIAA/NASA/ISSMO  Symposium  on  Multidisciplinary  Analysis  and 
Optimization,  Raphael  T.  Haftka  and  Layne  T.  Watson. 

•  Alumni  Award  for  Excellence  in  Research  (VPI&SU),  Layne  T.  Watson. 
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